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ABSTRACT 


This  thesis  considers  the  interoperability  of  recent  modeling  efforts  that  apply  constrained 
optimization  (combined  with  representations  of  system  function  and  management)  to 
assess  and  improve  the  operational  resilience  of  critical  infrastructure  (Cl)  systems  to 
disruptive  events.  We  implement  these  mathematical  models  using  the  Pyomo 
optimization  package,  which  is  built  on  top  of  the  Python  programming  language.  This 
computational  environment  provides  advantages  for  data  preprocessing  and  post¬ 
processing,  including  convenient  and  efficient  methods  for  manipulating  Cl  network  data. 
Moreover,  the  object-oriented  nature  of  Pyomo  creates  a  natural  means  for  representing 
interdependent  Cl  systems.  Specifically,  the  model  for  each  Cl  system  can  be 
implemented  as  its  own  object,  and  the  combined  model  can  be  implemented  as  another 
object  built  from  its  dependent  components.  This  allows  for  increased  flexibility  and 
extensibility  beyond  previous  implementations.  We  manage  the  inputs  and  outputs  of  the 
models  in  a  way  to  be  able  to  compare  them  across  studies,  obtaining  insight  on  their 
performance,  interactions,  and  effectiveness.  This  thesis  supports  a  broader  effort  to  build 
a  repository  of  functional  Cl  models  enabled  from  a  geospatial  user  interface  and 
connected  to  a  common,  backend  simulation  engine. 
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EXECUTIVE  SUMMARY 


Earthquakes,  floods,  typhoons,  and  climate  change,  as  well  as  terrorist  attacks, 
can  cause — and  have  caused — damage  to  critical  infrastructural  (Cl)  systems. 
Interruptions  of  energy  (electricity,  fuel,  gas)  supply  networks  or  communication  (glass- 
fiber-optic)  networks  can  be  damaging  because  of  our  dependence  on  them.  Moreover,  a 
disruption  in  one  Cl  system  can  have  a  cascading  effect  on  another  Cl  system  because 
these  systems  are  commonly  interdependent.  Because  of  the  desire  for  more  resilient  Cl 
systems,  there  is  an  ongoing  need  to  develop  and  apply  operational  models  of 
infrastructure  function. 

Unfortunately,  previous  research  and  development  (R&D)  efforts  to  improve 
infrastructure  resilience  have  been  hampered  because  of  a  lack  of  model  interoperability, 
a  lack  of  functional  realism,  sensitivity  of  real  system  geography  or  vulnerability,  and/or 
an  inability  to  visualize  system  behavior  in  the  presence  of  disruptive  events. 

This  thesis  addresses  the  interoperability  of  recent  modeling  efforts  at  the  Naval 
Postgraduate  School  (NPS)  that  apply  constrained  optimization  (combined  with 
representations  of  system  function  and  management)  to  assess  and  improve  the 
operational  resilience  of  Cl  systems  to  disruptive  events.  We  implement  these 
mathematical  models  using  the  Pyomo  optimization  package,  which  is  built  on  top  of  the 
Python  programming  language.  This  computational  environment  provides  advantages  for 
data  preprocessing  and  post-processing,  including  convenient  and  efficient  methods  for 
manipulating  and  visualizing  Cl  network  data.  Moreover,  the  object-oriented  nature  of 
Pyomo  creates  a  natural  means  for  representing  interdependent  Cl  systems.  Specifically, 
the  model  for  each  Cl  system  can  be  implemented  as  its  own  object,  and  the  combined 
model  can  be  implemented  as  another  object  built  from  its  dependent  components.  This 
allows  for  increased  flexibility  and  extensibility  beyond  previous  implementations. 

We  consider  the  behavior  of  a  notional  pair  of  interdependent  fuel  and  electricity 
infrastructure  systems.  We  manage  the  inputs  and  outputs  of  the  models  in  a  way  to  study 
them  in  isolation  or  as  a  combined  pair,  obtaining  insight  on  their  perfonnance, 


xv 


interactions,  and  effectiveness.  We  validate  our  results  against  those  from  a  previous 
study. 

This  thesis  supports  a  broader  effort  to  build  a  repository  of  functional  Cl  models 
enabled  from  a  geospatial  user  interface  connected  to  a  common,  backend  simulation 
engine. 
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I.  INTRODUCTION 


Earthquakes,  floods,  typhoons,  and  climate  change,  as  well  as  terrorist  attacks, 
can  and  have  caused  damage  to  critical  infrastructure  (Cl)  systems.  Our  complex  society 
highly  relies  on  these  systems,  and  therefore  an  effort  to  make  them  more  resilient  to 
natural  or  human  incidents  is  necessary.  Interruptions  of  energy  (electricity,  fuel,  gas) 
supply  networks  or  communication  (glass-fiber-optic)  networks  are  additionally  highly 
dependent  on  each  other.  Therefore,  optimizing  usage  and  structure  of  one  system 
realistically  means  optimizing  the  dependent  systems  appropriately  in  context.  Studies 
dealing  with  optimizing  networks  of  infrastructure  used  software  to  simulate  destructive 
events  and  improve  resilience. 

The  common  concept  of  operational  resilience  for  a  Cl  system  to  adapt  and  to 
maintain  its  function  in  the  presence  of  a  disruptive  event  is  introduced  in  “National 
Strategy  for  Homeland  Security”  (Homeland  Security  Council  [HSC]  2007),  which  was 
driven  by  the  Executive  Order  13010  establishing  the  the  U.S.  President’s  Commission 
on  Critical  Infrastructure  Protection  (PCCIP)  in  1996.  The  process  led  to  the  Presidential 
Policy  Directive  21  (PPD)  (White  House  2013),  and  the  concept  was  expanded  in  the 
National  Infrastructure  Protection  Plan  (NIPP)  in  December  2013  by  the  Department  of 
Homeland  Security  (DHS) 

Alderson  et  al.  (2014)  describe  how  to  build  and  solve  a  sequence  of  models  in 
order  to  assess  and  improve  the  operational  resilience  of  an  infrastructure  system  for 
worst-case  disruptions.  However,  in  order  to  support  real  decision-makers,  these  analyses 
require  two  important  features: 

•  Functional  Realism:  The  behavior  of  the  model  has  to  be  as  close  as 
possible  to  the  function  of  a  real  infrastructure  system.  For  example,  a 
model  of  electric  power  transmission  should  reflect  the  laws  of 
electromagnetics  (e.g.,  Kirchhoff  s  Laws). 
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•  Geographical  Realism:  Cl  systems  are  inherently  geospatial.  When 
superimposed  on  a  map,  they  should  look  as  close  as  possible  to  a  real 
system. 

In  practice,  however,  the  use  of  real  infrastructure  system  data  and  geography  for 
research  is  often  restricted  to  avoid  revealing  vulnerabilities  to  a  potential  attacker.  As  a 
result,  at  the  moment  there  is  an  inability  to  share  realistic  infrastructure  models  and  data, 
and  this  is  hampering  research  and  development  (R&D)  efforts  to  improve  infrastructure 
resilience.  Each  study  so  far  is  a  one-time  result,  which  does  not  allow  teaching,  results  in 
lost  insights  after  execution  and  presentation,  and  sometimes  leads  to  repetition  of  work. 
Additionally,  there  are  no  benchmarks  for  comparison  of  standard  data  sets  as  model 
inputs,  standard  functional  models,  canonical  examples  of  resilience  or  brittleness,  and 
algorithm  perfonnance. 

Overcoming  these  barriers  requires  that  we  address  several  questions:  How  can 
the  inputs  and  outputs  of  a  Cl  model  be  standardized  without  loss  of  generality?  Is  it 
possible  to  use  the  result  from  one  model  in  another  to  improve  the  results?  Is  there  a  way 
to  incorporate  different  models  so  that  solving  them  optimizes  both  models?  Is  there  a 
way  to  compare  two  models  that  analyze  the  same  system  and  get  a  measure  of 
effectiveness  in  a  way  that  hints  at  improving  them? 

This  thesis  considers  the  interoperability  analysis  of  recent  simulations  that  apply 
constrained  optimization  (combined  with  models  of  system  function  and  management)  to 
assess  and  improve  the  operational  resilience  of  critical  infrastructure  systems  to 
disruptive  events.  A  goal  of  this  research  is  to  study  interactions  between  models  and 
systems  and  to  measure  effectiveness  of  models. 

This  thesis  focuses  on  combining  the  results  of  resilience  analysis  of  critical 
infrastructure  using  network  models  via  a  set  of  software  solutions  that  integrate  different 
types  of  models.  We  manage  the  inputs  and  outputs  of  the  models  to  compare  them 
across  studies.  We  implement  the  results  from  one  model  into  another  using  backend 
simulation  engines  to  run  optimization  on  a  full-spectrum  set  of  influences.  The  goal  is  to 
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build  a  repository  of  functional  Cl  models  enabled  from  a  geospatial  user  interface  and 
connected  to  a  backend  simulation  engine  (see  Figure  1). 


Geospatial 

Data 


Economic  Data  Engineering  Data 


User  Interface 


Figure  1.  Integration  of  real-world  behavior  rules  and  operational 
infrastructure  models  via  a  single  computational  simulation  engine 
plus  visualization  and  interaction  with  it  via  a  GUI  (from  Naval 
Postgraduate  School  Center  for  Infrastructure  Defense  [CID]  2013) 


We  apply  attacker-defender  models,  as  in  Alderson  et  al.  (2014),  to  problems 
involving  interconnected  infrastructure  systems  in  the  spirit  of  Dixon  (2011)  and 
following  the  mathematics  developed  in  Dickenson  (2014).  We  consider  a  geospatial 
computational  platfonn  as  in  Martin  (2014).  The  results,  performance,  and  usability  are 
compared  to  the  onetime  results  from  previous  studies. 

The  ultimate  goal  of  the  research  is  to  create  a  sandbox  for  analysis,  design,  and 
virtual  testing  of  new  infrastructure  systems  and  operation  that  allow  what-if  scenarios  to 
system  responses  to  disruptive  events  (e.g.,  cascading  failure  behavior)  in  a  context-rich, 
controlled  environment  for  exercises,  education,  and  training  on  realistic  Cl  systems.  The 
hope  is  that  this  sandbox  serves  as  a  better  environment  to  produce  predictions  about 
unknown  interactions  of  the  grids  and  to  assess  resilience  for  real  systems. 
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Chapter  II  of  this  thesis  provides  a  brief  summary  of  previous  work  that  has  been 
done  and  shows  the  possibilities  for  improvement.  Chapter  III  presents  the  mathematical 
models  we  use  in  our  analysis.  Chapter  IV  shows  how  we  integrate  interacting  models 
into  a  structured  and  improved  computational  environment  that  enabled  the  user  to 
modify  their  tests  and  get  their  transferable  results  via  a  GUI  (graphical  user  interface) 
from  a  backend  server.  Chapter  V  summarizes  our  contribution  and  outlines  the  potential 
for  future  work. 
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II.  LITERATURE  REVIEW 


A.  PREVIOUS  WORK 

There  already  exist  several  models  and  algorithms  to  improve  networks  of 
infrastructure  in  terms  of  resilience.  The  literature  on  this  is  large  and  growing.  A 
tremendous  amount  of  work  has  been  done  by  the  Center  for  Infrastructure  Defense 
(CID)  at  the  Naval  Postgraduate  School  (NPS).  This  chapter  provides  an  abbreviated 
review  of  the  work  most  relevant  to  this  thesis,  along  with  an  outline  of  our  contribution. 

The  general  attacker-defender  (AD)  and  defender-attacker-defender  (DAD) 
models  for  Cl  are  defined  in  Brown  et  al.  (2006).  These  models  help  identify  functional 
and  important  parts  of  the  network  and  seek  a  solution  to  gain  resilience  by  using  game 
theory  and  optimization  models.  They  consider  deliberate  disruptions  of  functional 
networks  by  intelligent  adversaries.  They  assess  worst-case  effects  of  multiple 
disruptions  to  identify  the  system’s  ability  to  return  to  a  normal  state. 

Salmeron  et  al.  (2004,  2009)  apply  AD  models  to  the  electric  power  system.  They 
implement  dominant  physical  rules,  like  Kirchhoff  s  Laws,  in  a  functional  network  that 
therefore  allows  realistic  predictions  of  this  nontrivial  construct. 

AD  and  DAD  models  are  used  to  study  fuel  infrastructure  in  Alderson  et  al. 
(2015).  For  a  notional  Cl  system,  they  use  a  sequence  of  models  to  assess  and  improve 
operational  resilience,  which  then  can  be  achieved  by  actions  like  installing  redundant 
components  or  capacity  expansion. 

The  basic  mathematics  and  modeling  for  interdependent  infrastructure  systems 
were  created  by  Dixon  (2011).  Dixon  first  classifies  relationships  and  interdependencies 
in  the  analysis  of  effects  on  Cl;  he  then  uses  this  data  in  a  game-theory  AD  model  to  gain 
insight  into  otherwise  unpredictable  effects  caused  by  these  interactions.  Finally,  he 
provides  mechanisms  and  procedures  to  use  this  on  real-world  data  and  shows  how  to 
give  a  decision-maker  information  that  he  can  understand  and  use  with  network  analysis 
procedures  to  get  low-cost  suggestions  with  a  significant  improvement  in  overall  cost. 
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Dickenson  (2014)  expanded  the  work  of  Dixon  (2011)  to  a  larger  application 
modeling  interdependent  fuel  and  electric  power  infrastructures.  Dickenson  fonnulates  an 
algorithm  to  analyze  dependent  electric  and  fuel  infrastructure  systems  using  penalty 
functions  to  weight  the  interactions.  He  solves  for  the  combined  behavior  of  separate 
dependent  models  and  compares  the  results  to  those  from  a  single  integrated  monolithic 
model  and  finally  gives  suggestions  of  further  improvements  possible. 

Martin  (2014)  presents  a  geospatial  tool  for  infrastructure  analysis.  Martin 
develops  a  graphical  interface  to  make  it  easier  to  input  parameters  for  a  minimum  cost 
flow  network  problem.  He  focuses  on  the  analysis  of  a  fictional  Cl  system,  which  has  a 
geospatially  and  functionally  realistic  structure.  He  uses  this  tool  to  evaluate  the 
resilience  of  this  fictional  system. 

This  thesis  follows  the  basic  steps  for  assessing  and  improving  infrastructure 
resilience,  as  defined  by  Alderson  et  al.  (2014).  Their  tutorial  shows  how  to  build  a 
sequence  of  models  to  improve  the  resilience  of  CI.  Disruptions  can  be  caused  by  nature 
or  human  action.  They  use  simplified  examples  and  real-world  analogies  to  clarify  the 
impact  and  the  advantages  of  existing  models,  to  improve  the  state  of  being  prepared,  and 
to  enhance  new  structures. 

B.  OUR  CONTRIBUTION  IN  CONTEXT 

This  thesis  focuses  on  optimizing  the  results  of  resilience  analysis  of  critical 
infrastructure  using  network  models.  In  addition,  it  enables  interoperability  of  the  models 
and  gives  insight  into  the  performance  and  interactions  of  given  models.  The  goal  is  to 
build  a  repository  of  functional  CI  models  enabled  from  a  geospatial  user  interface  and 
connected  to  a  backend  simulation  engine. 
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III.  MATHEMATICAL  MODELS 


This  thesis  adopts  the  fuel  and  electric  networks  from  Alderson  et  al.  (2014)  and 
Sahneron  et  al.  (2004),  respectively,  to  formulate  a  combined  mathematical  model  as  in 
Dickenson  (2014)  to  allow  these  separate  networks  to  interact  in  tenns  of  suppliers  of  one 
system  demanding  a  commodity  of  the  other  system. 

A  key  building  block  in  our  analysis  is  a  mathematical  construct  used  by  Dixon 
(2011).  He  connects  supply  at  a  node  in  one  network  to  the  flow  on  an  arc  in  another 
network.  As  illustrated  in  Figure  2,  flow  K  on  arc  (/,/)  of  infrastructure  77  depends  on  a 

flow  Vnij  from  node  n  of  infrastructure  network  7  reaching  a  specific  threshold. 


Figure  2.  The  requirement  of  a  commodity  from  infrastructure  I  by 
infrastructure  II  (from  Dixon  2011). 


Mathematically,  this  dependence  is  represented  using  the  following  two 
constraints: 


threshold {j  ■  72  <  Vnjj 


(1) 


Y  <u.-T 

V  lJ  V 


(2) 


72  is  a  binary  variable.  It  is  used  to  determine  whether  the  threshold  requirement 
is  met.  72  is  set  to  zero  when  the  flow  Vmj  of  the  required  commodity  I  is  below  the 
requirement  threshold ^ .  72  can  be  set  to  one  only  when  Vnij  is  greater  than  or  equal  to 
threshold (j .  72  is  the  flow  in  infrastructure  77.  If  72  is  zero,  the  capacity  u;>.  of  arc  (/,/)  is 
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set  to  zero  and  no  flow  Yy  is  possible.  If  Ty  is  one,  then  the  capacity  ujy  of  arc  (ij)  is  the 


regular  upper  bound  of  the  flow  Yy . 

This  setup  is  used  on  all  dependencies  between  the  electric  and  the  fuel  network. 
The  networks  from  Dickenson  (2014)  are  not  altered  in  our  analysis  to  compare  results 
and  derive  insights  from  the  analysis. 

A.  SETS,  DATA,  AND  VARIABLES  FOR  OUR  FORMULATIONS 

For  completeness  and  consistency,  we  present  all  index  sets,  data,  and  variables 
for  our  formulations,  as  originally  described  in  Dickenson  (2014). 

1.  Indices  and  Index  Sets 


n  g  FN 

nodes  in  fuel  network  (alias 

n  g  FD  c=  FN 

demand  nodes  in  fuel  network 

neFS^FN 

supply  nodes  in  fuel  network 

FArcs  c=  FN  x  FN 

arcs  in  fuel  network 

n  g  PN 

nodes  in  power  network  (alias  pi,  pj) 

n  e  PD  c=  PN 

demand  nodes  in  power  network 

n  g  PG  <z  PN 

generation  (supply)  nodes  in  power  network 

n  g  PI  c=  PN 

bus  nodes  in  power  network  (where  supply  =  0) 

PArcs  c:  PN  x  PN 

arcs  in  power  network 

n  g  PDN  c  PN 

power  demand  nodes  that  supply  fuel  components 

(fi,  Jj)  e  PDFAn  c  FArcs 

power-dependent  fuel  arcs:  ( fi,jj )  g  PDFAn  can  carry  flow 

if  and  only  if  power  supply  to  n  g  PDN  exceeds  a  given 

threshold. 

n  g  FDV  c  F/V 

fuel  demand  nodes  that  supply  power  components 

[pi,  pj)  e  FDPAn  ci  PArcs  fuel-dependent  power  arcs:  ( pi,  pj)  e  FDPAn  can  carry 

flow  if  and  only  if  fuel  supply  to  n  e  FDN  exceeds  a 
given  threshold 
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2. 


Parameters  [Units] 


F Supply  n 
FLdSCostn 
FArcCost  fl  g 
FArcCap  fi  fj 
PDemn 
PArcCap  pi  pJ 
PYpi.pj 
PThreshn 

PGenCapn 

PGenCostn 

PLdSCostn 

PArcRespipj 
PArcReapipj 
B  . 

pi,pj 

FYfl.jp 

FThreshn 


fuel  supply  at  node  n  e  FN  [Barrel(bbl)/  hour(hr)] 
fuel  load  shedding  cost  of  demand  node  n  e  FD  [$/bbl] 
per-unit  cost  to  move  fuel  on  arc  ( fi,fj )  [$/bbl] 
capacity  of  fuel  on  arc  ( fi,fj )  [bbl/hr] 
power  demand  at  node  n  [Megawatts(MW)] 
power  capacity  of  ( pi,  pj )  [MW] 
power  flow  across  power  arc  ( pi,pj )  [MW] 

power  threshold  required  by  power  demand  node  n  e  PDN 
[MW] 

power  generator  capacity  of  n  e  PG  [MW] 

power  generator  cost  per  MW  by  node  n  e  PG  [$/Megawatt- 

hours  (MWh)] 

power  load  shedding  cost  of  node  n  e  PD  [$/MWh] 

resistance  of  arc  ( pi,pj )  [ohms] 

reactance  of  arc  (pi,pj)  [ohms] 

susceptance  of  arc  (pi,pj)  [l/ohms] 

fixed  fuel  flow  across  arc  ( fi,fj )  [bbl/hr] 

fuel  threshold  required  by  power  generation  node  n  e  FDN 
[bbl/hr] 


3.  Decision  Variables  [Units] 

FYfif]  flow  on  fuel  arc  (  /],  //  ]  [bbl/hr] 

FLdSn  load  shedding  at  fuel  demand  node  n  e  FD  [bbl/hr] 
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FTn 

PGenn 
PY  . 

p‘.pj 

PLdSn 

On 

PTn 


f  1  if  net  supply  to  fuel  demand  node  n  e  PDN  meets  or  exceeds  threshold 
[O  otherwise 

power  generated  at  generator  node  n  e  PG  [MW] 

flow  on  power  arc  ( pi,pj)  [MW] 

load  shedding  at  power  demand  node  n  e  PD  [MW] 

phase  angle  at  power  node  n  e  PN  [radians] 

f  1  if  net  supply  to  power  demand  node  n  e  FDN  meets  or  exceeds  threshold 
0  otherwise 


B.  FUEL  MODEL  WITH  INTER-DEPENDENCE 

We  begin  with  a  formulation  for  a  fuel  infrastructure  system,  whose  pumps 
depend  on  electric  power  provided  by  a  separate  infrastructure. 


2.  Discussion  on  Inter-dependence 

Key  to  the  functional  inter-dependence  are  constraints  (F3),  (F4),  and  (F6) 

involving  the  variables  FTn .  Each  FTn  is  a  switch  modeling  the  interdependence.  It  is  set 

to  one  if  the  net  electric  supply  received  by  demand  node  n  within  the  electrical 
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distribution  system  meets  or  exceeds  the  power  requirement  or  threshold,  PThreshn . 
Constraint  (F3)  sets  the  capacity  of  an  arc  of  the  fuel  network  to  zero  or  full  capacity,  if 
the  dependence  threshold  variables  FTn  are  set  to  zero  or  one.  Constraint  (F4)  sets  the 

dependence  threshold  variable  FTn  based  on  the  operating  conditions  in  the  electric 
system.  In  this  model  the  FTn  variables  are  only  constrained  by  the  fixed  electric  flows, 
PYpj  n ,  from  that  system.  Stipulation  (F6)  sets  the  dependence  threshold  variables,  FTn ,  as 
binary. 


C.  ELECTRIC  MODEL  WITH  INTER-DEPENDENCE 

We  continue  with  a  formulation  for  an  electricity  infrastructure,  whose  generators 
depend  on  the  energy  resources  of  a  separate  fuel  infrastructure. 


1.  Formulation 

min  Y  PGenCost nPGenn  +  V  PLdSCostnPLdSn 

PLdi 

s.t. 


PY,PGen  “G 
PLdS 


nePD 


pY,i,,i=Bpi.Pi(0pi'en) 

£  £  pr,.»=° 


pi:(pi,n) 

EPArcs 


Pj(n,pj) 

EPArcs 


pGen„+  X  PYrL1-  £  PY,,0=O 


pi:(pi,n) 
e  PA  res 


pj-{n,pj) 

EPArcs 


£  PYr,,~  £  PY,,r,+ppdSn=Pdem„ 


y  (pi,  pj)  e  PArcs 
Vn  e  PI 

Vn  e  PG 

Vn  e  FF 


pi:(pi,n) 

EPArcs 


Pj-(n,pj) 

EPArcs 


FThreshnPTn<  £  FY/,>  -  £  FY„,fj  VneFDN 


eFDPN 


eFDPN 


PYpi,pj  <  PA  rcCap pi  pj  PTn 
0  ^  pYpiPj  <  PArcCap pipj 
0  <  PGenn  <  PGenCapn 
0  <  PLdSn  <  PDemn 

pt.^M 


F n  e  FDN,  ( pi,  pj )  e  FDPAn 
y  (pi, pj)  e  PArcs 
V«  e  PG 
yn  e 

Vn  e  FD/V 


(PO) 


(PI) 

(P2) 

(P3) 

(P4) 

(P5) 

(P«) 

(P7) 

(PS) 

(P9) 

(P10) 
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2.  Discussion  on  Inter-dependence 

Key  to  the  functional  inter-dependence  are  constraints  (P5)  and  (P6)  involving  the 
variables  PTn .  Each  PTn  is  a  switch  modeling  the  interdependence.  It  is  set  to  one,  if  the 
net  fuel  supplied  to  demand  node  n  in  the  fuel  distribution  system  meets  or  exceeds  its 
fuel  requirement  (threshold),  FThreshn .  Constraint  (P5)  sets  the  dependence  threshold 

variable  PTn  based  on  the  operating  conditions  in  the  fuel  model.  In  the  fuel  model  the 

PTn  variables  are  only  constrained  by  the  fixed  fuel  flows,  FY ,  from  that  system. 
Constraint  (P6)  sets  the  capacity  of  an  arc  of  the  electric  network  to  zero  or  full  capacity, 
if  the  dependence  threshold  variables  PTn  are  set  to  zero  or  one. 

D.  COMBINED  MODEL 

We  combine  the  formulations  of  the  electricity  infrastructure  and  the  fuel 
infrastructure  into  a  single,  integrative  model,  and  convert  the  flow  parameters  in  the 
dependency  constraints  to  variables  in  order  to  solve  for  optimal  solution  of  the  combined 
model. 


1.  Formulation 


min  I  FA rc Cos t fi  /;  F Yfi  tj  +  ^  FLdSCostnFLdSn 

FY  ,FLdS  (fi,Jj)eFArcs  heFD 

FT,PY,PY 

PGen,PLdS 

+  ^  PGenCost nPGenn  +  y  PLdSCostnPLdSn  (CO) 

nePG  nePD 

s.t. 

(F1),(F2),(F3),(F5),(F6) 

(P1),(P2),(P3),(P4),(P6),(P7),(P8),(P9),(P10) 


PThreshnFTn  <  £  PYp„ 

-  y  py  . 

»,PJ 

V  n e PDN 

(Cl) 

pi:(pi.n ) 

Pj\n,pj ) 

ePArcs 

e^PArcs 

FThreshnPTn  <  £  FYf, 

i 

M 

Vn e FDN 

(C2) 

eFArcs 

fr-(n,S) 

eFAvcs 
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2.  Discussion  on  Inter-dependence 

Constraints  (Cl)  and  (C2)  model  the  interdependencies  between  the  two  systems. 
Constraints  (Cl)  requires  that  net  supply  of  electrical  power  to  power  demand  nodes  meet 
or  exceed  the  threshold,  PThreshn ,  to  allow  fuel  flow  on  the  corresponding 
interdependent  fuel  arc  Constraints  (C2)  requires  that  net  supply  of  fuel  to  fuel 

demand  nodes  meet  or  exceed  the  threshold  requirement,  FThreshn ,  to  allow  electric 
flow  on  the  corresponding  interdependent  power  arc  ( pi,pj ) . 

E.  DISCUSSION 

Dickenson  (2014)  implemented  this  model  using  the  General  Algebraic  Modeling 
System  (GAMS  2015),  but  his  implementation  consists  of  many  separate  model  files  (see 
Figure  3)  and  is  not  easy  to  extend  or  modify.  The  results  from  that  implementation  are 
simple  text  files  with  no  fixed  format,  which  makes  them  hard  to  use  in  common 
statistical  analysis  tools  like  JMP  (2015)  or  R  (2015) 


Figure  3.  GAMS  (.gms)  code  files  necessary  to  solve  the  combined  model 
and  the  corresponding  input  text  and  comma  separated  value  (.csv) 
files  as  implemented  in  Dickenson  (2014). 
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Building  on  Dickenson’s  basic  implementation,  we  make  several  improvements  to 
allow  extension  and  modification  of  the  model  and  allow  implementation  of  different 
solvers  and  more  advanced  algorithms. 
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IV.  IMPLEMENTATION  AND  RESULTS 


A.  OUR  APPROACH 

Dickenson  (2014)  demonstrates  how  two  well-understood  infrastructure  networks 
with  established,  published,  and  proven  formulations  could  interact  in  a  GAMS  (2015) 
environment.  The  first  one  is  taken  from  Alderson  et  al.  (2014),  a  fuel  model  with  their 
operator  formulation.  The  other  is  the  IEEE  (Institute  of  Electrical  and  Electronics 
Engineers)  RTS-96  (Reliability  Test  System  1996)  with  the  electric  model  taken  from 
Sahneron  et  al.  (2004).  The  interdependence  is  modeled  using  the  constructs  from  Dixon 
(2011). 

We  elect  to  implement  these  mathematical  models  using  the  Pyomo  optimization 
package  (2015),  which  has  been  built  on  top  of  the  Python  programming  language 
(2015).  There  are  several  reasons  for  this.  First,  a  Python-based  programming 
environment  has  advantages  for  data  preprocessing,  and  the  availability  of  libraries,  such 
as  NetworkX  (2015),  provide  convenient  and  efficient  methods  for  manipulating  Cl 
network  data.  Second,  Python  supports  subroutines,  whereas  GAMS  does  not,  and 
subroutines  are  essential  for  modular  algorithm  development.  Third,  the  object-oriented 
nature  of  Pyomo  creates  a  natural  means  of  representing  Cl  systems.  Specifically,  the 
operator  model  for  each  Cl  system  can  be  implemented  as  its  own  object,  and  the 
combined  model  can  be  implemented  as  another  object  built  from  its  dependent 
components.  This  allows  for  increased  flexibility  and  extensibility  beyond  similar 
implementations  in  GAMS.  Finally,  a  Python  implementation  allows  for  easier 
integration  into  other  applications,  specifically  the  Graphical  User  Interface  (GUI) 
defined  in  Martin  (2014),  which  was  also  implemented  in  Python.  Specifically,  we  are 
able  to  connect  our  models  to  that  GUI  via  a  transfer  file  using  the  extensible  Markup 
Language  (XML)  fonnat  (see  Figure  4). 
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Figure  4.  GUI  improves  usability  as  the  translation  into  one  program 
language.  One  type  of  transfer  file  that  is  universally  readable 
improves  exchange  and  analysis  option  of  solution  sets. 


B.  FOCUS 

The  goal  is  to  allow  a  continuous  and  reusable  integration  of  mathematical  models 
and  geospatial  data  into  a  graphical  representation  supported  by  software  that  is 
extensible  and  low-cost.  This  way,  all  work  can  be  reused,  further  improvements  and 
deeper  insights  can  be  gained,  and  there  is  no  need  to  reinvent  a  system  or  to  relearn  or 
reprogram  already  established  and  proven  results  when  using  them  on  a  bigger  scale  or  in 
different  model  environments. 

C.  MILESTONES 

We  subdivided  the  project  into  three  milestones  (A,  B,  and  C;  see  Figure  5). 
Milestone  A  involves  using  the  logic  of  Dickenson  (2014)  for  the  separate  fuel  model 
(FM)  and  electric  model  (EM)  in  Python  Pyomo.  Milestone  B  involves  building  a 
combined  model  (model  EM-FM)  in  an  efficient  way  that  produces  the  same  results  as 
the  original  GAMS  code.  Milestone  C  involves  using  the  model  from  the  GUI  and 
improving  further  aspects  of  usability,  security,  transferability,  and  efficiency.  We 
present  results  from  model  FM,  model  EM,  and  model  FM+EM  in  three  representative 
cases  to  compare  resulting  solution  network  flow  and  objective. 
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Figure  5.  Milestones  A,  B,  and  C  and  its  phases  for  our  translation  and 
improvement  process. 

1.  Milestone  A 

Milestone  A  is  a  logical  model  in  a  Python  Pyomo  environment  that  implements 
all  necessary  logical  steps  from  the  FM  in  GAMS  and  the  EM  in  GAMS. 

We  create  a  Python  subroutine  that  uses  Pyomo  to  build  the  FM,  which  includes 
all  sets,  all  parameters,  all  variables,  all  constraints  (including  dependency  constraints), 
and  the  objective,  as  defined  in  Section  III.A  of  this  thesis. 

We  build  on  the  notional  fuel  infrastructure,  which  is  introduced  in  Alderson  et  al. 
(2015)  and  shown  in  Figure  6  as  basis  for  our  analysis. 
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Figure  6.  Fuel  distribution  network  model  where  the  black  filled  circles  are 
the  supply  nodes  and  the  others  have  fuel  demands;  the  arcs 
represent  establish  connections  between  nodes  (from  Alderson  et 
al.  2015). 


We  create  another  Python  subroutine  that  uses  Pyomo  to  build  the  EM,  which 
includes  all  sets,  all  parameters,  all  variables,  all  constraints  (including  dependency 
constraints),  and  the  objective,  as  defined  in  Section  III.  A  of  this  thesis. 

Figure  7  shows  the  electric  system.  It  has  74  total  nodes  composed  of  33 
generator  nodes,  17  demand  nodes,  and  24  bus  nodes  (from  Reliability  Test  System 
[RTS]  Task  Force  1999). 
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Figure  7.  The  IEEE  One  Area  RTS-96  is  the  basis  for  the  electrical  power 
system  Salmeron  et  al.  (2004). 

2.  Milestone  B 

Milestone  B  is  a  combined  Model  (FM+EM)  that  consists  of  the  two  logical 
models  and  the  real  data.  It  allows  us  to  solve  the  two  models  separately  and  in  a 
combined  manner. 

We  create  a  Python  subroutine  that  uses  Pyomo  to  build  the  combined  logical 
model  (FM+EM),  which  imports  logical  models  FM  and  EM,  defines  links  between 
models,  and  returns  a  separate  logical  model.  Figure  8  shows  the  steps  and  inclusions  per 
model. 
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♦♦additional  constraintsto  link  EM  to  FM 

♦♦combined  objective  -^logical  model  FM+EM 

Figure  8.  Logical  models  of  fuel  model,  electric  model,  and  combined  model 
and  their  inclusions. 


We  create  separate  instances  for  model  FM,  model  EM,  and  combined  model 
FM+EM.  This  allows  us  to  solve  FM  alone  and  to  solve  EM  alone,  but  also  to  solve 
FM+EM  independent  of  the  single  models. 

a.  Comparison  of  Example  Used  in  Dickenson  (2014) 

The  interaction  options  are  shown  in  Figure  9  (electric  network  depends  on  fuel) 
and  Figure  10  (fuel  network  depends  on  electricity)  and  are  used  in  the  same  manner  in 
the  analysis  of  the  interdependency  study  from  Dickenson  (2014). 
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Figure  9.  RTS-96  positions  of  fuel  dependence  (Dickenson  2014)  where  the 
generator  nodes  demand  fuel  supply.  The  dependency  is  depicted 
in  dashed  lines.  For  example,  electric  flow  on  arc  (gl  1 5,i  1 15) 
depends  on  fuel  from  node  fill. 


Figure  10.  The  Fuel  Demand  Model  positions  of  electrical  power  dependency 
(Dickenson  2014)  where  pumps  depend  on  electricity.  The 
dependency  is  depicted  in  dashed  lines.  For  example,  fuel  flow  on 
arc  (fnl0,fn6)  depends  on  electricity  from  node  dl03. 
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The  next  two  figures  sum  up  the  output  and  the  results  from  the  run  of  the 
combined  model  FM+EM  in  GAMS  (see  Figure  11)  and  in  Python  Pyomo  (see  Figure 
12).  They  show  the  state  of  the  objective  values  for  the  separate  models  after  completion 
of  the  combined  run  in  the  structure  used  by  Dickenson  (2014).  The  baseline  solutions  of 
the  three  models  are  EM:  4460.3,  FM:  79593.5,  and  FM+EM:  84053.8.  These  models  use 
the  same  input  files  for  all  implementations.  The  specific  values  of  flows  differ  between 
the  implementations,  but  that  is  caused  by  the  fact  that  there  are  multiple  optimal 
solutions  for  the  flows. 
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Figure  1 1 .  Output  of  flows  and  results  from  the  base  case  for  model  FM+EM 
as  implemented  in  GAMS. 
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Figure  12.  Output  of  flows  and  results  from  the  base  case  for  model  FM+EM 
as  implemented  in  Python  Pyomo. 

b.  Comparison  of  Two  Examples  of  Stress  on  the  Inter-dependency 

We  simulate  three  situations  to  test  the  models  and  verify  their  results.  The  first 
one  checks  the  functionality  of  the  implementations,  and  the  following  two  examples 
show,  first,  a  case  in  which  the  GAMS  and  the  Python  Pyomo  version  have  to  find  a 
single  best  solution  to  reroute  the  commodity  of  one  network  with  high  cost  to  satisfy 
demands. 
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(1)  Case  Zero:  No  Load  Shedding  Cost 

We  set  the  load  shedding  cost  for  both  models  to  zero,  as  indicated  in  Table  1.  We 
expect  a  reduction  in  total  cost.  We  are  testing  basic  functionality  of  both  models  with  a 
known  outcome. 

Table  1 .  Load  shedding  costs  in  base  case  and  in  the  test  model. 


Model 

FM 

EM 

FM+EM 

load  shedding  cost 
before 

150 

150 

150 

load  shedding  cost 
after 

0 

0 

0 

The  effect  is  what  we  would  expect  because  the  optimization  for  minimal  cost 
decides  that  there  is  no  need  to  transport  or  satisfy  demands  at  all.  This  was  a  critical 
finding  in  Dickenson  (2014)  because  the  search  for  an  appropriate  high  penalty  for  load 
shedding  was  crucial  to  gaining  realistic  results.  These  results  match  the  two 
implementations  and  concur  with  the  expected  outcomes  (see  Figure  13). 


Python 
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Objective  from  EM: 

0 

Objective  from  FM+EM: 
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GAMS 

Objective  from  FM: 

0 

Objective  from  EM: 

0 

Objective  from  FM+EM: 

0 

Figure  13.  Results  of  case  zero  in  GAMS  and  Python  models. 


(2)  Case  One:  Electric  Generator  Out  of  Order 

In  the  first  example,  we  assume  an  electric  generator  (gl  18)  has  malfunctioned 
and  produces  zero  MWh  instead  of  400,  as  depicted  in  Table  2.  We  expect  no  influence 
on  the  FM,  but  a  cost  increase  in  the  EM  and  the  FM+EM,  because  the  total  power 
demand  cannot  be  satisfied  and  load  shedding  cost  should  occur.  We  are  testing  the 
optimization  process  of  both  systems  with  an  artificial  test  and  a  known,  but  non-ob vious 
solution. 
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Table  2.  Electric  generator  gl  18  has  capacity  of  0  instead  of  400  MWh. 


Model 

FM 

EM 

FM+EM 

capacity 

before 

- 

400 

400 

capacity 

after 

- 

0 

0 

All  three  models  produce  the  same  optimal  results,  and  the  effect  is  only  in  the 
EM  and  the  FM+EM  models.  This  test  shows  that  both  implementations  solve  the  single 
and  the  combined  model  independent  from  each  other.  The  penalty  in  EM  is  equal  to  the 
penalty  in  EM+FM  because  the  system  is  still  solvable  in  the  way  to  avoid  an  additional 
penalty  in  the  interdependent  FM  via  the  optimal  routing  in  the  separate  EM,  which  was 
artificially  chosen  to  test  the  optimization  process.  These  results  also  match  the  two 
implementations  and  concur  with  the  expected  outcome  tendencies  (see  Figure  14). 
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Figure  14.  Electric  generator  is  not  producing  energy  and  causes  a  17% 
increase  in  cost  in  the  EM. 


(3)  Case  Two:  Fuel  Arc  Broken 

In  the  second  case,  we  assume  a  fuel  arc  (FN8-FN12)  is  broken  and  has  zero 
capacity  instead  of  1350  bbl,  as  depicted  in  Table  3.  We  expect  an  influence  on  the  FM 
because  the  highest  supplier  FN8  has  to  reroute  its  fuel  to  satisfy  the  demands  within  the 
network.  We  do  not  expect  an  effect  on  the  interdependent  electric  network  because  the 
amount  of  fuel  that  can  be  transported  has  not  changed.  We  are  testing  the  optimization 
process  of  both  systems. 
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Table  3.  Fuel  arc  FN8-FN12  capacity  of  0  instead  of  1350  bbl. 


Model 

FM 

EM 

FM+EM 

capacity 

before 

1350 

- 

1350 

capacity 

after 

0 

- 

0 

All  three  models  produce  the  same  optimal  results,  and  the  effect  is  only  in  FM 
and  FM+EM.  Every  arc  in  FM  has  the  same  capacity  and  is  set  so  that  the  highest 
supplier  is  able  to  send  his  product  completely  over  any  arc.  This  test  shows  that  even  in 
this  comfortable  case,  the  rerouting  cost  to  avoid  penalty  forms  demands  within  the 
system  and  penalty  from  the  interdependent  system  can  get  relatively  high  (see  Figure 
15). 
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Figure  15.  Broken  fuel  arc  causes  a  26%  increase  in  cost  in  FM. 


(4)  Discussion 

Both  GAMS  and  Python  Pyomo  models  give  the  same  objective  values  and  are 
able  to  solve  the  problems  in  a  short  amount  of  time.  An  evaluation  of  computation  time 
can  be  performed  only  with  more  systems  to  be  compared  and  with  systems  with  more 
nodes  or  a  higher  complexity.  Extreme  cases  occurring  in  fuel  and  in  electric  network  on 
dependent  nodes  cause  results  represented  in  both  models  as  intuitively  suspected  or 
artificially  designed  to  test  their  reliability  in  function. 

3.  Milestone  C 

Milestone  C  consists  of  several  improvements  that  seemed  necessary  in  terms  of 
usability,  security,  transferability,  and  usability  of  a  backend  simulation  engine. 
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We  adapted  the  GUI  originally  developed  from  the  Infrastructure  Design  Editor 
Analyzer  (IDEA)  project  (Johnson  and  Heine  2015)  to  work  in  a  Microsoft  (2015b). 
Windows  8.1  and  therefore  future  systems.  We  managed  the  usage  of  a  single  backend 
server  system  to  allow  restriction  of  access  to  sensitive  data. 

We  redraw  the  fuel  network  Alderson  (2014)  into  the  IDEA  (Johnson  and  Heine 
2015;  see  Figure  16)  environment  and  programmed  Python  Pyomo  (2015)  code.  This 
model  uses  a  CPLEX  (2015)  solver  but  can  technically  use  different  solvers  by  changing 
only  a  single  line  of  code. 

An  additional  adaption  was  to  store  variations  of  the  solution  sets  in  a  reasonable 
manner  that,  for  a  specific  network,  the  user  is  able  to  switch  between  different  input 
settings  and  immediately  see  the  visualization  on  screen  while  also  being  able  to  extract 
the  data  for  use  in  an  arbitrary  tool  for  further  analysis  or  representation  (JMP  2015, 
Microsoft  2015a,  R  2015).  The  networks  can  be  created  and  modified  in  the  IDEA  Editor 
(see  Figure  16). 


Figure  16.  Screenshot  of  the  IDEA  Editor  showing  the  creation  of  the  basic 
fuel  network. 


The  interaction  between  these  two  networks  can  now  be  executed  with  a  single 


click  on  a  “Run  Simulation”  button  after  loading  them  into  the  GUI.  The  flows  are  color 
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coded  and  can  also  be  varied  by  the  user,  as  depicted  in  Figure  17.  The  results  proved  to 
be  equal  to  those  obtained  in  Dickenson  (2014). 


Figure  17.  Screenshot  of  the  IDEA  Analyzers  showing  a  solution  for  the  basic 
fuel  network  without  background  in  preset  color  coding  for  flows 
and  the  ability  to  compare  results  visually  and  via  tables  quickly. 


In  integration  of  a  Dystopia  (Center  for  Homeland  Defense  and  Security  [CHDS] 
2011)  network,  layers  from  the  provided  server  worked  out  well,  and  an  integration  of 
real-life  data  is  possible.  Figures  18  and  19  show  the  artificial  island  of  Dystopia  with 
rudimentary  networks  drawn  already. 
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Figure  18.  The  Dystopia  map  from  the  Center  for  Homeland  Defense  and 
Security  (2015)  and  its  developed  networks  so  far. 


Figure  19.  Screenshot  of  a  Dystopia  map  used  by  Martin  (2014)  to  show  a 
partial  network  to  visualize  results  to  the  background.  The  example 
here  is  a  telecommunications  network  but  is  illustrative  of  the 
ultimate  goal  to  connect  models  of  interdependent  Cl  systems  in  a 
geographically  realistic  environment. 


To  make  integration  in  the  user  interface  easier,  we  included  a  help  file  and  a 
step-by-step  example  to  the  IDEA  Editor  so  that  the  usage  of  this  tool  behaves  more  like 
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any  other  standard  GUI  and  allows  quick  access  to  commonly  asked  questions  and  a 
convenient  introduction  for  first-time  users.  The  structure  is  shown  in  Figure  20  (without 
the  27  steps  of  the  tutorial). 
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Figure  20.  The  structure  of  the  IDEA  Editor  help  file  allows  users  to  take  a 
tutorial,  find  their  topic  on  one  site,  and  use  the  search  option 
effectively. 
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V.  SUMMARY  AND  CONCLUSION 


A.  SUMMARY 

Our  goal  has  been  to  create  a  user-friendly,  adaptable  basis  for  future  analysis  on 
interdependent  Cl  networks.  We  accomplish  several  steps  in  that  process.  We  replace 
multiple  programming  language  solutions  (GAMS,  MS  Excel,  etc.)  with  a  single 
freeware  language  (Python  Pyomo).  In  the  model,  we  allow  any  solver  to  be  used  and 
show  the  output  with  CPLEX.  It  is  now  possible  to  work  from  anywhere  via  a  server  on 
any  recent  Microsoft  Windows  system.  The  distribution  of  information  is  now 
controllable  due  to  the  backend  solution.  The  intermediate  files  are  reduced  from  over 
twenty  to  two  and  standardized  in  one  globally  usable  fonnat  (XML).  We  code  a 
functional  Cl  fuel  model.  We  code  a  functional  Cl  electric  model.  We  combine  these 
models  in  a  new  model.  Editing  and  visualisation  via  a  GUI  are  now  possible.  Analysis 
becomes  faster,  more  user-friendly,  and  in  accordance  with  notional  standards.  Analysis 
modifications  in  the  GUI  are  also  more  user-friendly,  and  their  outcomes  are  provided  in 
transferable  files.  The  final  models  can  be  used  to  improve  resilience  of  networks  and 
find  vulnerable  spots  in  case  of  a  possible  disaster.  This  way,  a  blackout  or  loss  of 
communication  becomes  less  likely,  and  the  chance  of  losing  a  high  amount  of  human 
life  can  be  reduced.  Overall,  our  contribution  allows  easier  access  and  faster  analysis  of 
interdependent  Cl  networks. 

B.  CONCLUSION 

We  transfer  two  non-trivial  logical  network  models  and  their  data  into  a  single 
open-source  programming  language  to  allow  future  work  to  be  implemented  more  easily. 
We  make  results  transferable  and  enable  the  user  to  interact  with  the  models  via  a  GUI 
that  gives  visual  feedback  on  results  and  has  the  opportunity  to  change  settings  and 
values  quickly.  With  all  the  other  improvements  stated,  we  found  a  base  for  future  studies 
that  are  not  one-time  results  anymore,  and  we  enhanced  attractiveness  for  using  it  in 
education,  research,  and  presentations. 
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c. 


FUTURE  WORK 


To  improve  global  usage,  there  is  a  need  to  develop  models  for  other  networks 
that  might  interact  (e.g.,  telecommunication  or  transportation)  with  realistic  behavior  and 
a  clear  dependency  link  between  the  systems. 

To  improve  the  usage  of  our  tool,  we  should  expose  students  of  OR,  researching 
faculty,  managers  of  Cl,  and  decision-makers  of  cities/industry/military  to  the  system  and 
find  out  whether  the  visualization  can  be  improved.  We  should  also  create  a  stack  of 
more  models  and  networks  to  do  research  and  out-of-the-box  experiments  with  the 
networks  that  are  not  real  but  represent  the  behavior  of  each  system  and  its  interaction  in 
a  realistic  way. 

To  improve  in  tenns  of  user-friendliness,  it  would  be  useful  to  create  a  tool 
similar  to  the  editor  and  the  analyzer  that  allows  for  creating  scenarios  without  the  need 
to  manually  enter  every  change  into  the  data  or  the  GUI,  but  rather  set  up  rules  to 
influence  the  network  in  cases  of  disruptive  events  like  an  earthquake,  tsunami,  or  bomb 
attack,  for  example,  via  a  separate  Scenario  Editor  GUI  (see  Figure  21). 
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Figure  2 1 .  Future  situation  of  data  organization  and  applications  including  the 
Scenario  Builder  and  a  library  of  verified  GIS  (Geographic 
Information  System)  data. 
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Because  our  system  allows  integration  of  different  algorithms  in  a  convenient 
manner,  we  recommend  testing  and  exploring  different  approaches,  which  is  optimal  for 
large-scale  networks  and  a  higher  number  of  interacting  networks.  For  example,  to 
validate  the  results,  a  remote  and  known  network  of  a  feasible-sized  infrastructure,  like 
the  electric  and  fuel  power  grid  of  the  Hawaiian  Islands,  could  be  used  in  future  studies. 
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